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Abstract We have developed a mathematical model for in-host virus dynamics that 
includes spatial chemotaxis and diffusion across a two dimensional surface repre- 
senting the vaginal or rectal epithelium at primary HIV infection. A linear stability 
analysis of the steady state solutions identified conditions for Turing instability pat- 
tern formation. We have solved the model equations numerically using parameter 
values obtained from previous experimental results for HIV infections. Simulations 
of the model for this surface show hot spots of infection. Understanding this local- 
ization is an important step in the ability to correctly model early HIV infection. 
These spatial variations also have implications for the development and effectiveness 
of microbicides against HIV 
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1 Introduction 

The normal response in an individual after being infected by a virus is activation of 
the immune system, driving infection levels down. If the immune response is suffi- 
ciently potent then the disease can be completely eradicated from the body, but in 
many instances this does not occur. Instead, over the course of time, an eventual bal- 
ance of disease replication and immune clearance is established leading to chronic 
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infection. These steady state outcomes, clearance versus chronic infection, are sug- 
gestive of simple dynamics but this ignores spatial variations, including possible hot 
spots of infection (5], which confound the dynamics in both transient and long time 
behaviours. Although there has been some mathematical modelling of acute HIV in- 
fection I2II I2 II28L current mathematical models of HIV infection mainly focus on 
response to antiretroviral therapy of HIV viral levels after the viral setpoint 12711311 . 
Moreover these models assume a well-mixed environment with no real spatial be- 
haviour. This is very different to what happens at the very earliest stages of infection 
at the vaginal or rectal epithelium during sexual transmission of HIV from an infected 
man to his partner. A single lineage usually expands in the new host, even though the 
genetically heterogeneous inoculum contains numerous infectious units. The dynam- 
ics of high HIV seminal loads leading to sporadic infection and the establishment of 
single foci of infection 11111 . are difficult to understand biologically and completely 
fall outside the sphere of usual mathematical modelling of infectious diseases with 
simple ordinary differential equations. 



Spatially heterogeneous outcomes may arise from underlying spatial heterogene- 
ity possibly due to tissue architecture or damage from other sexually transmitted in- 
fections [4 1, but it is also possible to have such non-uniformities spontaneously arise 
from the infection dynamics. One of the archetypical manners in which this occurs 
is through the existence of a Turing instability in a set of partial differential equa- 
tions (PDEs). Turing's seminal work in 1952 [29 1 showed that for some nonlinear 
reaction-diffusion equations the steady state solution of the system is not spatially 
uniform. The Turing instability occurs when a spatially homogeneous steady state of 
the reaction dynamics, which is linearly stable in the absence of diffusion, becomes 
linearly unstable when the reactions are coupled with the diffusion. This can occur 
when there are two or more nonlinearly interacting species with different diffusivi- 
ties. The resultant inhomogeneous spatial pattern is called a Turing pattern. Turing 
patterns have been proposed to explain patterning in numerous physical, chemical 
and biological systems (TJ, including models of morphogenesis [19,3, IT31[T7l and 
some chemical reactions [25]. Turing pattern formation has also been investigated in 
an SIR model to predict the spatial transmission of diseases in a population [ 14 1. 



In order to investigate the impact of spatial dynamics in a simple mathemati- 
cal model of HIV infection we extended an SIR model for in-host virus dynamics 
B24II23I to include spatially random diffusion and spatially directed chemotaxis 1101 . 
The spatio-temporal behaviour of this system is investigated within the framework of 
Turing instability-induced pattern formation [29] and is shown to result in spatial hot 
spots of infection for certain ranges of the parameter values. We further explored the 
behaviour of the model system through numerical simulations which reveals compli- 
cated spatial dynamics persisting through transient and long time behaviours. These 
spatial variations have implications for the development and effectiveness of micro- 
bicides against HIV (H). 
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2 Model Equations 

The standard SIR-based model for in-host virus dynamics is given by 131U23ll27ll 

dT 

— =s-kVT-uT 
dt 

^-=kVT-8l (1) 
dt 

dV 

— =NI-cV. 
dt 

In these equations the dependent variables are: T, the population of uninfected target 
cells; I, the population of infected cells; and V, the population of free virions, where 
these all vary with time t, but not space. It is assumed that the target cells are supplied 
at a constant rate s and they are removed either through cell death with death rate jj, 
or by becoming infected by virions. The parameter k represents the rate of infection 
of target cells per virion; N represents the number of virions produced per unit time, 
per infected cell; 8 is the death rate of the infected cells and c is the clearance rate of 
virions. 

Extending this model to also include spatial aspects so that the independent vari- 
ables are now (f,x) produces the following model 

d T 

— =s - kVT - fiT + D T V 2 T - #V(rV7) 
at 

d ^=kVT-8l + D I V 2 I (2) 
dt 

— =NI-cV+D v V 2 V, 
at 

where T, I and V are now concentrations of target cells, infected cells and virus 
respectively, with appropriate units according to the space dimension (e.g. cells/mm 
for ID or cells/mm 2 for 2D). In this model it is assumed that the target CD4+ T cells, 
infected cells and free virions all diffuse with diffusion constants Dt, D[ and Dy re- 
spectively. We have also included a spatial chemotaxis term — %V{TVI) to represent 
the chemotactic attraction of target immune cells driven by the concentration gradient 
of cytokines from inflammation at sites of infection. The random walk diffusive mo- 
tion of these T cells has been well established both in vitro and in vivo 1 18 1 with some 
evidence for an anomalous component to the diffusion [6 1, although in this work, for 
simplicity, we will assume a purely Brownian diffusion. The chemotaxis of T cells 
is more difficult to establish but careful experiments have clearly demonstrated the 
chemotaxis of T cells in response to gradients of chemokines in microfluidic in vitro 
studies lfT3ll . The vaginal or rectal epithelium will be represented by a 2-dimensional 
surface so that our space component is given by x = (xi,X2) T . 



3 Turing patterns 



We are interested in whether the dynamics of the spatially extended model allow for 
formation of spatial patterns. The system of partial differential equations |2]) may 
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be classified as a reaction-chemotaxis-diffusion system. Patterns may occur in this 
system in the neighbourhood of a spatially homogeneous steady state provided the 
conditions for a Turing instability are met, namely that this spatially homogeneous 
steady state is: 

(Tl) linearly stable in the absence of diffusion and chemotaxis; and 
(T2) linearly unstable in the presence of diffusion and chemotaxis. 



4 Steady States 

The spatially extended system Q, has the same (spatially) homogeneous steady 
states as the standard model ([TJ: the disease-free state 

r * = -i jg = o, v *=o, 

M 

and the endemic state 

r * = c5 skN-cSn = skN-c8fi 

W' Nk8 ' ck8 



5 Non-dimensional equations 



In order to reduce the number of parameters and simplify some of the analysis it 
is useful to work with a non-dimensional version of the system |2]). Let the non- 
dimensional dependent variables be given by 

u\ = T/T c 

U 2 = I/I c 
«3 - V/V c 

and the non-dimensional independent space and time variables be given by 

Xi =xi/L 
X 2 =x 2 /L 

T = t/t c 

where the values of T C ,I C , V C ,L and t c shall be chosen later in such a way to minimise 
the total number of parameters. We substitute the non-dimensional variables into ([2]) 
and obtain 



^ = (~) - (v c t c k)u 3U1 - (fit c )m + v 2 «i - v ( Ml v« 2 ) 

du 2 (kV c T c t c \ {WArt 
-j— = I — I uiu 3 - (dt c )u 2 + I I v "2 

du 3 ( NI c t c \ ( D v t c \ 2 

^ = V^r) U2 ~ { ctc)U3+ {-ir) VM3 ' 
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There is now freedom in setting the scaling factors and simplifying the equations. 
To this end set T c , I c and V c to the corresponding values of the endemic homogeneous 
steady state T*,I* ,V*. The non-dimensional equations then become 

du\ ( skN\ . . / Djt c \ _2 

(1 -U3U1) -t c flm(l ~U3)+ I — y- V U 



. — t c 1 . 1 \ ■ n.iii 1 1 1,7.1 11 j v • 11 , 1 1 , , 

-<«(< f ^)™ 

^ u 2 , c \ , <j s f Dit c \ _ 2 

= (dt^uxiiT, - (dt c )u 2 + I I V mi 



(4) 



<9m3 

^7 



(cf e )M 2 - (cf c ) M 3 + (^-JY^j 



V 2 M 3 . 



We may further simplify |4} by choosing f e := and L 2 := / fi. Then we set 
the new non-dimensional parameters to be the following 

skN D l D v 

<,:--?—, dr.— —, dy := — , 

s% ( 1 A 5 c 



Finally, we arrive at the non-dimensional version of Q: 

-^r- - = ^ — — l)u\Uj — u\ + V 2 mi — c/y V (mi Vm 2 ) 
<9m2 . . 9 

— — = a(MiM3 — m 2 ) +rf/V m 2 (j) 

5m3 2 
-= — = p (M2 — M3 ) + dy V M3 . 

It is straightforward to check that in the absence of spatial variation this system admits 
the following two equilibria: 

(El) The endemic spatially homogeneous steady-state at (m^m^Mj) = (1,1,1); and 
(E2) The disease-free spatially homogeneous steady state at (m*,^,^) = (1,0,0). 

Observe also that the value of ^ in |5| determines one of two distinct scenarios: 

(i) When | < 1 we have I,V < whenever m 2 ,M3 > 0, hence the endemic homoge- 
neous steady state is not physically relevant as it corresponds to negative concen- 
trations of infected cells and free virus. The only spatially homogeneous steady 
state is the state free of disease, (^ ,0,0). 

(ii) When £, > 1 both of the spatially homogeneous steady states, (^,0,0) and (1,1,1), 
are physically relevant. 

In order to explore the possibility of Turing pattern formation we shall consider 
the stability of |5) linearised about each of the two steady states. 
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6 Linearising about the endemic spatially homogeneous steady state 

Let us assume for this section that ^ > 1 as we have already established above that 
this is a necessary condition for a physically relevant endemic homogeneous steady 
state. We perturb this equilibrium (u^u^u^) = (1,1,1) by writing u\ = u* + Aui for 
each z € {1,2,3}. The linearised form of |5]) for the perturbations Alii is given by 



d_ 



Au\ 








Au\ 




"1 -d x " 




Au\ 


Au2 




a -a a 




Au2 


+ 


d l 


V 2 


Au2 


Aui 




_0 j8 -p 




Auj, 




d v 




Aut, 



(6) 



It shall be convenient to carry out Fourier transforms with respect to the spatial vari- 
ables in equation (|6j. This yields 



d_ 



~AU{ 




AU 2 




AU 3 





\-i;-q 2 ) d x q 2 



a 




{-a-diq 2 



-(1-1) 
a 



dyq 2 ) 





~AU{ 




AU 2 




AU 3 



(7) 



where AUj denote the Fourier transform of Auj and q 2 := q T q where q £ R d is the 
Fourier variable, with d the spatial dimension. The conditions for Turing instabilities 
are determined from the eigenvalue spectrum of the coefficient matrix in |7]). The 
requirement that the homogeneous steady state is stable in the absence of diffusion 
and chemotaxis is met if all eigenvalues have negative real parts when q 2 = 0. A 
Turing instability may then occur if one or more eigenvalues have positive real parts 
for some q 2 > 0. 

The characteristic polynomial of the matrix in (|7| is 

?, 3 + X 2 (E,+a + p + (l+d I + d v )q 2 )) 
+X {£ > a + E > fi + {ad v +fid I + E,d I + a + E,d v +fi-ad x )q 2 
+ {did v +di + d v )q 4 } (8) 
+ { (4 ad v + £,j5d] - aj5d x )q 2 + (£,d]d v + ad v + fidi - ad x dy)q 4 
+ (did v )q 6 + {^-\)ap} 

To check for stability, recall that a cubic polynomial A 3 + ak 2 + bX + c has all 
roots in the left half complex planeif and only if a, b, c > and ab> c (Ruth-Hurwitz 
stability criterion for a cubic polynomial |9|). Below we consider the characteristic 
polynomial for different values of q 2 . 



6.1 Case of no spatial variation 

In the absence of diffusion and chemotaxis (q 2 = 0) the characteristic polynomial ([8]) 
simplifies to 



A 3 +A 2 (a + /3+£) + A(c^+j3£) + (£-l)a/3. 



(9) 
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Since we have assumed that S, > 1, all the coefficients of (|9]l are positive. Thus 
the remaining condition 

(a + j3 + £)(a + j3)<^ >(|-l)a/3 (10) 

is necessary and sufficient for all roots of ^ to be in the left half-plane. 
Since a and j3 are positive and S, > 1, equation (JTOj is always true: 

>(S-l)ap. 

Thus the endemic steady state is stable in the absence of diffusion and chemotaxis 
and the first Turing condition (Tl) is satisfied. 



6.2 Case of spatial variation 

Now we investigate the nature of the roots of (|5]l in the presence of diffusion and 
chemotaxis (q 2 > 0). 

Firstly, consider the case for large q 2 . Then the characteristic polynomial <(8j will 
have all its coefficients positive (as in each term the largest power of q 2 has a positive 
coefficient). Moreover, the coefficient of A 2 tends to (1 + dj + dy)q 2 , the coefficient 
of A tends towards (d/dy + dj + dy)q 4 while the constant coefficient is dominated 
by (djdy)q 6 . It is then straightforward to show that the product of the former two is 
larger than the latter, hence for sufficiently large q 2 all roots of ([8]) are in the left half 
complex plane and the system is stable. 

Even though for large values of q 2 the system becomes stable, there may still be 
sufficiently large values of d x and an appropriate range of values of q 2 for which 
instabilities occur. To determine the threshold on d x (or %) above which instabilities 
may be possible, let us first write the polynomial ([8]) in the following form: 

A 3 + A 2 («! + a 2 q 2 ) + A(&i + b 2 q 2 + b 3 q 4 ) + (ci + c 2 q 2 + c 3 q 4 + c 4 q 6 ) (11) 

for appropriate values of coefficients ai,a2,bi,b2,b^,ci,C2,c 3 and C4. Note that all 
these coefficients are positive except for possibly b 2 , c 2 and C3. If we assume that 
these too are positive then: 

(51) aibi > c\ 

(52) aib 2 +a 2 bi > c 2 

(53) aibj +a 2 b 2 > C3 

(54) a 2 bi > C4. 
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The proof of this is straightforward but tedious, and is shown in the Appendix. Now, 
(S1)-(S4) imply that 

(a\+a 2 q ){b\+b 2 q + b 3 q )>c\+c 2 q + c 3 q +c 4 q , 

for all q, therefore by the Ruth-Hurwitz stability criterion, all roots of (jTTJ (and (|8]l) 
are in the left half complex plane and thus |7]i is stable. 

We conclude that the only way the system in Q may become unstable is if at least 
one of b 2 , c 2 or c 3 becomes non-positive. Therefore the necessary (but not sufficient) 
condition for Turing condition (T2) to hold is that either 

(CI) b 2 <0,ihati&d x >dy + (P + Z)d I /a + l + % + l/a; 
(C2) c 2 < 0, that is d % > t;d v /J3 + %di /a; or 
(C3) c 3 < 0, that is d % > \di/a + 1 + f}d[/(ad v ). 

Using d x = (l - ^ we see that (CI), (C2) and (C3) are equivalent to 

(CI') X > q( f gfa (dv + 03 + g )di/a + 1 + g + l/a), 

(C2') X> sx ^) (W<i + Wa),or 

(C3') X > , z ff [ /{) ^//« + 1 + Pdi/(ad v ) 

respectively. Given a set of parameters, for the purpose of determining a possible 
existence of a Turing instability we shall only consider the weakest condition of the 
three. 



7 Linearising about the disease-free spatially homogeneous steady state 

Similar to our approach for the endemic equilibrium in the previous section, we may 
linearise |5]) about the disease-free spatially homogeneous steady state (u^,u^,u^) = 
(<i;,0,0). Since this state is always physically relevant, for the time being we need 
not assume any additional conditions on | (apart from positivity). After performing 
a Fourier transform of |5]) linearized about (§,0,0), we get 



~AUi~ 




~-l-q 2 d z $q 2 


-1(1-1)" 




~AUi~ 


AU 2 




(-a-diq 2 ) 


at 




AU 2 


AU 3 




j8 


(-V-dyq 2 ) 




AU 3 



In the spatially homogeneous setting (q = 0), it is easy to check that the disease- 
free steady state is stable if and only if £, < 1 . 

The matrix in (jT2j has — 1 — q 2 < as one of its eigenvalues. The remaining two 
are given by the eigenvalues of the 2x2 lower-right submatrix 

(-a-diq 2 ) at; 

P (-p-d v q 2 ). 

It is easy to check that if § < 1 the eigenvalues of this matrix lie in the left half 
complex plane (e.g. observe that the matrix has negative trace and positive determi- 
nant), and this holds for all values of q 2 . So (Tl) and (T2) cannot both be true, and 
we conclude that the Turing conditions cannot be satisfied in a neighbourhood of the 
disease-free steady state. 
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8 Regularisation of chemotaxis 

It is well-known that in two spatial dimensions chemotaxis above a certain threshold 
results in a finite time blow up of solutions to the governing equations (see e.g. [8]). 
The standard approach to deal with this non-realistic phenomenon is to introduce a 
regularisation term to (|5j. We shall use a density-dependent sensitivity regularisation, 
studied in |30|, based on the assumption that with increasing cell density, their ad- 
vective velocity reduces. For other forms of regularisation see the survey article by 
Hillen and Painter [7|. Our new governing equations become: 

= £-(£ - 1)ki^-ki+V 2 mi -(1+eWyV ( — — — V M2 

dx \\+EU\ 

*^ = a{mu\ -u2) + diV 2 u2 (13) 
dx 

-zr^- = fi(ll2 — M3) +£A/V 2 M3 

OX 

for some dimensionless regularisation parameter e > such that in the limit as e — > 
we recover the original non-regularised model. 

Let us introduce the notion of effective chemotaxis: 

^ (Mi;e):= TW^ 

Then d x = d x at the endemic steady state (u\ = 1), and the linearised equations |6]) 
and Q remain unchanged, therefore in this case our stability analysis and conditions 
for Turing patterns formation from previous sections also apply in the regularised 
setting. Also note that d x — > as u\ — > °°. 



9 Parameters, units and dimensions 

There has been considerable investigation of parameter values for variants of the 
SIR model, based on trials in HIV-infected individuals H27ll26ll20ll22l and from SIV 
infected macaques lfl6l . and it is reasonable to assume that these parameter values 
provide useful starting approximations for variants of the model Q, that also includes 
a spatial component. 

We assume the original HIV model ([T]) has parameter estimates given by jV = 
480 virions cell -1 day -1 , k = 3.43 x 10 ml virions -1 day -1 , 8 — 0.5 day -1 , c = 
3 day -1 , s = 10 cells mm -3 day -1 , /J. — 0.03 day -1 . 

We will perform numerical simulations in both one and two spatial dimensions. 
Note that the constants k and s given above are volume-based. In order to adapt them 
to two or one spatial dimensions let us assume that the region in which we are solv- 
ing the PDE is either a thin rectangular sheet of thickness h or a thin wire with a 
square hxh cross-section. Thus the new dimension-specific values of k and s become 
k = k\r'^ d and s — s/z 3-rf respectively, where d G {1,2} is the number of spatial di- 
mensions. Note that by changing d the only non-dimensional parameter that changes 
is d y . For all of the numerical simulations we shall take h = 0.1 mm. 
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The diffusion of T cells in lymphatic tissue has been estimated at 

D T = 1.1 nm 2 s -1 =0.09504 mm 2 day" 1 (14) 

It can be assumed that the uninfected {T) and infected (7) CD4+ T cells will have 
similar diffusion coefficients. Different studies have shown virions with diffusion co- 
efficients of the order of 0.0088jtim 2 s . The only parameters without experimental 
bounds are the effective chemotaxis term % and the regularisation constant e. For 
the above stated parameters, the weakest necessary (but not sufficient) condition for 
Turing instability is (C2), requiring d x > 219.8, that is % > 10.4 mm 4 cell -1 day -1 
in two spatial dimensions and % > 104 mm 3 cell -1 day -1 in one spatial dimension. 
These are our lower bounds on the chemotactic threshold. 

Figure [T] shows the real part of the leading eigenvalue of the matrix of the lin- 
earized system |7]i for the parameters stated above as a function of spatial frequency. 
Note that even though ^ = 110 mm 3 cell -1 day -1 is above our calculated lower bound 
on the chemotactic threshold, the system remains stable for all spatial frequencies. 
Turing conditions are only met for % somewhere between 110 and 120 mm 3 cell -1 day 

Recall that we set the scaling factors for independent variables to t c = 1 //x and 
L 2 = Dj / ji. With the given parameters, this results in values t c « 33.3 days and 
Lw 1.78 mm. 



10 Numerical solutions in one spatial dimension 

The model PDE is solved numerically on a one-dimensional domain of length L with 
zero-flux boundary conditions. In one-dimension the chemotaxis term need not be 
regularised, thus, unless otherwise stated, we shall use e = in this section. We used 
MATLAB's built in ID PDE solver pdepe. The initial conditions are taken as either 
a small random or deterministic perturbation around the endemic steady state. 

For x below 104 mm 3 cell -1 day -1 (our lower bound on the chemotactic thresh- 
old for a Turing instability) the solutions become flat with no spatial features as time 
progresses. This can be seen in Figure|2] 

For x above this threshold, we see Turing pattern formation in the form of three 
uniformly separated peaks. This pattern remains and the peaks do not blow up in 
time; see Figure [3] 

As the strength of chemotaxis increases, we see that the spikes become steeper. 
However the solutions are still finite for all time; see Figure|4] 

Increasing the length of the domain results in proportionally more spikes; see 
Figure [5] It should also be noted that in this case it takes a longer length of time for 
the solutions to settle to a steady state and the spikes to become uniform in height. 

Starting with a random initial condition of ui (e.g. uniform on (0.975, 1.025)), 
the oscillations are initially at a higher frequency than for a smooth initial condition. 
However these oscillations quickly settle and produce a pattern similar to the case 
with normal initial distribution; see Figure [4] 
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Fig. 1 Leading eigenvalues of the linear model in one spatial dimension as a function of spatial fre- 
quency for several values of %. Empty squares correspond to % = 100mm 3 cell -1 day - '; empty cir- 
cles to % — 110mm 3 cell -1 day -1 ; empty triangles to % = 120mm 3 cell -1 day -1 ; filled squares to 
X = 130 mm 3 cell -1 day -1 ; filled circles to x = 150 mm 3 cell -1 day -1 ; and filled triangles to % = 
170 mm 3 cell -1 day -1 . The upper right panel shows an enlarged view around the origin. 



10.1 Dependence on initial conditions 

The spatial frequency of the steady state solution may depend on the initial conditions 
for long domains. Figure [6] shows density profiles for two solutions with different 
initial distributions of infected cells — the first one concentrated on the left border 
and second one concentrated in the middle of the region. The resulting final states 
contain a different number of peaks — 16 and 17, respectively. 



10.2 Solutions with initial conditions near the disease-free steady state. 

We have already shown that the disease-free homogeneous steady-state is unstable 
when skN > cdjJ.. Hence we can expect any disturbance of this state (i.e. the intro- 
duction of virus or infected cells) to result in system transitioning to the endemic 
steady state. 

Because we are considering global behaviour of the system it is more useful to 
draw plots in terms of re-dimensionalised variables (for example, it enables one to 
compare the number of infected cells to the number of healthy cells). In the absence 



12 



O. Stancevic et al. 
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Fig. 2 Solution for % = 100 mm 3 cell -1 day -1 . The solid blue line represents ui (left axis), the non- 
dimensional version of concentration of healthy cells and the dashed red line is uj (right axis), the con- 
centration of infected cells. The concentration of infected cells wi is initially normally distributed. The 
concentration of free virus (not shown in the graph) is initially constant at 143 = 1 and changes similarly to 

u 2 . 



of spatial variation a typical transition to the endemic steady state is shown in Figure 

m 

Figure[8]shows the transition from a small random disturbance around the disease- 
free steady state to the endemic steady state with chemotaxis % = 130 mm 3 cell -1 day 
above the critical threshold, regularised with e = 1 .0. The resulting Turing pattern is 
identical to the pattern produced starting from a perturbation of the endemic steady 
state (see e.g. Figure [3}. 



10.3 Consequences of changing the virus clearance rate 



Figure 10 shows the change in the Turing pattern after the rate of virus clearance in- 
creases from c = 3 day -1 to c = 4 day -1 . Even though the overall number of infected 
cell reduces, the maximum density of infected cells remains roughly the same as with 
the lower rate of clearance. Note that increasing c in this way will result in a lower 
chemotactic threshold for Turing pattern formation. 
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Fig. 3 Solution for % = 1 30 mm 3 cell - 1 day - 1 with same initial conditions as those in Figure]^ The solid 
blue line represents u\ (left axis), the non-dimensional version of concentration of healthy cells and the 
dashed red line is ui (right axis), the concentration of infected cells. 



11 Numerical Solutions in two spatial dimensions 



The two-dimensional model exhibits a blow up of solutions in finite time for suf- 
ficiently strong chemotactic attractions. For this reason we have employed a regu- 
larisation as described in (JT3J. Setting the thickness of the surface on which we are 
modelling the infection to h = 0.1mm, results in the critical value for chemotaxis % 
to lie somewhere between 1 1 and 12 mm 4 cell -1 day -1 . 

The governing PDE was solved numerically on a square domain with sides of 
length L « 1.78 mm (i.e. X\ = X2 = 1). The method of lines was used, which in- 
volved semi-discretising ( fTJI in both space variables on a uniform rectangular grid 
and solving the resulting ODE system in time using MATLAB's implementation of 
the Runge-Kutta method, ode45. The results for the endemic steady state are shown 



in Figure 1 1 and Figure 12 



While the two-dimensional model is more difficult to solve numerically and some 
form of chemotactic regularisation is required we did not observe any qualitative 
behaviour that was different from the setting of only one dimension in space. 



12 Discussion 



We have shown analytically and numerically that HIV infection need not be spatially 
homogeneous and can exhibit Turing patterns. The patterns can only occur if i) the 
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Fig. 4 Solution for £ = 220 mm 3 cell -1 day -1 with ui = «3 = 1 initially constant and «2 taken from a 
uniform random distribution in (0.975, 1.025). Note that the size of the domain has been doubled com- 
pared to the preceding figures. The solid blue line represents U\ (left axis), the non-dimensional version of 
concentration of healthy cells and the dashed red line is m (right axis), the concentration of infected cells. 



conditions for existence of an infected steady state hold (skN > c8fi), and ii) the 
chemotactic attraction is strong enough. With the parameter values as outlined in 
Section 9, the estimated time for transition to the nonhomogeneous endemic state to 
occur locally (i.e. on a small patch of tissue) is around 14 days (see Figure[7|. 

Larger chemotactic attraction results in a larger range of possible spatial frequen- 
cies for the Turing patterns as well as higher amplitudes of the peaks. While the 
resulting pattern may depend on the initial state of the infection, we found that in 
most cases it does not. However, different initial states may significantly affect the 
speed at which the pattern is established, with small random perturbations of the en- 
demic state enabling faster settling to a Turing pattern than localised perturbations 
which in turn settle faster than from a perturbation of the disease-free steady state. 

We found that if the initial perturbation is local, the pattern tends to get estab- 
lished near the perturbation, and then propagates outwards, whereas for random ini- 
tial perturbations the patterns emerge more or less simultaneously across the whole 
region. 

Upon getting infected by HIV it typically takes 1-2 months ETI for the body's 
immune system to respond by increasing the rate at which the virus is cleared. This 
is enough time for the initial Turing pattern to form. After the clearance rate is in- 
creased, it is possible for new, more prominent and less frequent Turing patterns to 
form. Also, even if the system was in a state that does not admit any patterns and the 
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is estimated to be 0.36 mm day -1 . The solid blue line represents u\ (left axis), the non-dimensional version 
of concentration of healthy cells and the dashed red line is U2 (right axis), the concentration of infected 
cells. 



infection becomes spatially homogeneous, such a response from the immune system 
may give rise to pattern formation. 

This analysis indicates that foci of HIV infection that are observed in tissue, can 
be established as a result of the dynamics of the system irrespective of any spatial 
heterogeneity in the tissue. These foci may provide an environment where new in- 
fection of cells outweighs immune clearance and hence supports the maintenance of 
infection despite an expanding immune response. Our analysis indicates that these 
patterns are established over a longer time scale than would be relevant for the im- 
pact of a microbicide which acts at the very earliest stages of infection. This analysis 
however assumes that the surface itself representing the vaginal or rectal epithelium, 
is homogeneous. This is not the case however and further work to assess the impact 
of the generation of Turing patterns in a more realistic environment is required. 
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Appendix: Showing conditions for Turing instability 

Proposition 1 Consider the characteristic polynomial \\\) . Provided a\,a2,b{,b2,b^,c\, C2,cj and C4 < 
all positive, conditions (S1)-(S4) hold- 
Proof We demonstrate each of the conditions in turn below. 
(SI) 

>t;ap 

>(§-l)a)J=ei; 



(S2) 



(S3) 



(S4) 



aib 2 + a 2 b\ =(% + a)b2+pb 2 + (l+di)bi+d v bi 
={£,+a)b 1 + (\+d l )b l 

+ fi ( ad v +Pdi + %di + a + £dy+p- ad x ) + i, ad v + £,fid v 
={£, + a)b 2 + (q + d,)bi +p(ad v + fid, + a + %d v + ft) + %pd v 

+ (^ad v + ^d I -aPd x ) 
>£,ad v +t;Pd] - afid x = c 2 ; 



a\b3+a 2 b 2 =(| + a + P)(d,d v +di+dv) 

+ {l + di + d v )(ad v + Pdi + %di + a + ^d v + P- ad x ) 
>t;d]d v + ad v +Pd] 

+ (1 + di)b 2 +dv(p2 + ad x ) - ad x d v 
>i s d[dy + ady + pdj — ad x dy = cy, 



a 2^3 =(1 +di + dy)(didy +di+dy) 
>d/dy = C4. 
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